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Recently, Wang and Landau proposed a new random walk algorithm that can be very efficiently 
applied to many problems. Subsequently, there has been numerous studies on the algorithm itself 
and many proposals for improvements were put forward. However, fundamental questions such as 
what determines the rate of convergence has not been answered. To understand the mechanism 
behind the Wang-Landau method, we did an error analysis and found that a steady state is reached 
where the fluctuations in the accumulated energy histogram saturate at values proportional to 
[log(/)]~^''^. This value is closely related to the error corrections to the Wang-Landau method. We 
also study the rate of convergence using different "tuning" parameters in the algorithm. 



I. INTRODUCTION 



Computational methods have been used extensively for solving complex problems in the past decades. In particular, 
in statistical physics equilibrium quantities of a system with many degrees of freedom are measured. The framework 
of statistical physics is formalized such that all equilibrium quantities can be derived from the partition function, 

Z(T)=^e-^('^)/'=-^ (1) 

a is the state of the system, E is the energy corresponding to cr, is the Boltzmann constant and T is the temperature. 
The summation is over all possible states and the number of possible states is a colossal number which cannot 
be enumerated. Nevertheless, computational methods such as Monte Carlo techniques P] are used to sample the 
partition function; in particular. Metropolis importance sampling 0] has achieved considerable success. However, 
new techniques are emerging and are replacing the Metropolis importance sampling especially near phase transition 
boundaries where the Metropolis importance sampling becomes inefficient. A class of new techniques, called the 
generalized ensemble methods, such as the multicanonical method 0,01, the broad histogram method 0] and the flat 
histogram method Q , were developed based on re- writing the partition function as a sum over energies 

Z{T) = e-^^"^/'"^ = [e-^('^)/'=-^] (2) 

{a} E 

where the partition function is reduced from a sum over all states to a sum over r^N energy levels. The partition 
function would be tractable if the energy density of states q{E)co\x\A be calculated. 

Recently, a systematic, iterative, random walk method H,|3^3 been proposed as one of the generalized ensemble 
methods. Now generally known as the Wang-Landau algorithm, it has received much attention in literature and has 
been applied to a wide range of problems 1 1. 12. 13. Mil. T here have also been numerous proposed improvements and 
studies of the efficiency of this method [TBI flall^llSl Il9l f20l l2ll 123. l2 ^ : however, there are still many unanswered 
questions, e.g. what determines the rate of convergence to the true density of states and is there any universality 
behavior related to this algorithm? In this paper, we attempt to quantify the mechanism behind the Wang-Landau 
method and study the effects of using different "tuning" parameters. 

The Wang-Landau algorithm ItA llOjl is an iterative process in which the density of states g{E) is modified by a 
factor /fc > 1, and the refinement of the density of states is assured with monotonically decreasing modification factors, 
e.g. fk+i = \fjk- For each time the energy level E is visited, g{E) is multiphed by /fc > 1, and a histogram of energy is 
accumulated concurrently. It was proposed that the modification factor be decreased when the accumulated histogram 
satisfies a certain flatness condition which we call the stopping condition. In this paper we study the effects of different 
modification factors and stopping conditions, and derive an expression for the error term in the Wang-Landau method 
based on generalizations of the modification factors and stopping conditions. We shall consider arbitrary sequences 
of modification factors, /i > • • • /fc > /fc+i • • • > I and arbitrary sequences of corresponding stopping conditions 
Ai, A2, • • •• The stopping conditions Ai, A2, • • • may or may not be the histogram flatness condition; other stopping 
conditions could be used, for example, stopping after some predetermined maximum number of Monte Carlo steps, 
or stopping after some number of times the random walker reaches the minimum energy state. This generalization is 
needed for theoretical error analysis in the next section. 
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II. THEORETICAL ERROR ANALYSIS 



In the Wang-Landau algorithm, for each visit to an energy level E, the density of states at that energy level increases 
by a factor fk > ^, and the corresponding histogram increases by one. Assume that initially the unknown density of 
states were set to 1, i.e. go{E) = IV E. The density of states at the end of the nth iteration is given by 

n 

logg„(i5) = ^i7,(i?)log(/fe) (3) 

k=l 

where Hk{E) is the accumulated histogram and fk is the modification factor for the fcth iteration. At this point, 
it is important to realize that the relative values of g{E) are sufficient for calculating thermodynamics quantities. 
Hence, a constant factor can be extracted from gn{E), without losing any information, by a change of variable on the 
histograms, 

Hk{E) ^ Hk{E) - xmxi{Hk{E)} ^ Hk{E) for fc = l,2,---,n (4) 

E 

where mhiE{Hk{E)^ is the minimum value of the accumulated histogram for the fcth iteration. Then Eq. Q becomes 

n 

log5„(£') = ^ Hk{E) log(/fe) + constant of energy. (5) 
fc=i 

The second term in Eq. |SJl is independent of energy, and from here onwards we shall refer to gn{E) as the density of 
states without the second term in Eq. (O, i.e. 

n 

\og~gn{E)^Y.Hk{E)\og{fk). (6) 
fc=i 

To lay the foundation for deriving an expression for the error of the Wang-Landau method, we use the conjecture 
that the method converges to the true density of states with proper choice of parameters. 

Conjecture 1 Let the Wang-Landau algorithm he carried out with a sequence of modification factors, ■ ■ ■ fk > fk+i > 
• • • /oo = 1- There exists a sequence of stopping conditions Ai, A2, • • • Aoo such that, 

lim gn{E) — goc{E) — g*{E) x constant (7) 

n — *oo " 

where gn{E) is the density of states calculated up to the nth iteration and g*{E) is the true density of states. 

This conjecture does not give any error bounds on the density of states; it only says that, in the limit of an infinite 
number of iterations the Wang-Landau estimate converges to the true density of states. In addition, no constraint is 
imposed on the stopping conditions in the conjecture. The error term up to the nth iteration can be defined as 

oc 

Y,[\og~gUE)-\ogUE)]^Y. E Hk{E)\og{fk) (8) 

E E k=n+l 

An intuitive view of Eq. (jHJ is that, if an infinite number of iterations were performed, the exact answer would be 
obtained. When n iterations were done instead, the error of gn{E) will be the sum of all the rest of the iterations 
that were not carried out explicitly. Define 

AHk = Y,HkiE) (9) 

E 

where AiJ^ is a measure of fluctuations in Hk{E). By the assumption of convergence series (implied by Conjecture 
the order of summations in Eq. ^ can be rearranged. Then, Eq. ^ becomes 

oc 

77„ = 5][log.9oo(i?)-logg„(^^)] = E Ai7fclog(/fe) (10) 

E k=n+l 

Eq. (|10|l shows that, assuming appropriate stopping conditions, 77„ depends only on the fluctuation in the histogram 
and the sequence of modification factors fk- When the values of fk are predetermined (e.g. fk+i = \/Jk); 
becomes the only determining factor for 7y„. 
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FIG. 1; A_fffc versus Monte Carlo steps per site for 16 x 16 ferromagnetic Ising model for various log(/) values. From left to 
right, top to bottom, log(/) values are 10"^, 10"^, 10"* and 10~^. 



III. RESULTS 



We investigate the Monte Carlo time dependence of Ai/fc for each iteration with the Wang-Landau method, where 
the subscript k denotes the fcth iteration. Simulations were performed on the ferromagnetic Ising model and on the 
fully frustrated Ising model with various system sizes. The Hamiltonian is 

n^-Y.Jv'^^'^j (11) 

where the sum is over nearest neighbors on a two dimensional square grid and ai takes the values ±1. Jij — 1 for 
the ferromagnetic Ising model, and for the fully frustrated Ising model, Jij takes the value —1 for every alternate 
horizontal nearest neighbors bonds and +1 otherwise. 

Fig. n shows the time dependence of AHk for several values of log(/); log(/) — 10"^, 10"^, 10"'' and 10~^ from 
left to right, top to bottom respectively. We used the sequence of correction factors log(/A:+i) = log(/fe)/1.78 with 
log(/i) = 0.1, this sequence is chosen so that log(/fc+4) = log(/A;)/10. These graphs were generated by performing the 
Wang-Landau algorithm on a 16 x 16 ferromagnetic Ising model with numerical values averaged over 128 independent 
simulations. The Monte Carlo steps per spin, the horizontal axis of Fig. ^ are measured from the time when we 
decrease log(/). AHk increases initially and eventually saturates, and for smaller log(/) values, saturation values 
are greater and number of Monte Carlo steps required to reach saturation are larger. Because the error term given 
by Eq. (|10|l depends only on AHk, any computational effort after AHk become saturated does not improve the 
accuracy of the final density of states gn{E). On the other hand, stopping the random walk before AHk becomes 
saturated would make the simulation less efficient because insufficient statistics are accumulated in the fcth iteration 
and much more statistics would have to be accumulated with smaller log(/) values for subsequent iterations. An 
optimal algorithm is to stop the simulation as soon as AHk becomes saturated. The Wang-Landau algorithm in the 
original paper suggested using the histogram ffatness condition as a stopping condition, but this does not guarantee 
optimal efficiency. 

It is difficult to predict the saturation value of AHk for fc = 1, • • •. As shown in Fig. ^ we performed a set of 
simulations with more Monte Carlo steps than required for AHk to reach saturation. In this way, we could measure 
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FIG. 2: Plots of saturation AHk versus log(/) for 8x8 ferromagnetic Ising model (filled circles), 8x8 fully frustrated Ising 
model (empty circles), 16 x 16 ferromagnetic Ising model (filled triangles) and 16 x 16 fully frustrated Ising model (empty 
triangles). 128 independent simulations were performed for each data point and error bars were smaller than the size of the 
symbols. 

the saturation values accurately. Fig. |21shows a plot of saturation values versus log(/) for ferromagnetic Ising model 
(FMIM) and fully frustrated Ising model (FFIM). In double log scale, the data points fall on a straight line with the 
values of the slopes equal to -0.491 ± 0.004 for 8 x 8 FMIM, -0.501 ± 0.004 for 8x8 FFIM, -0.496 ± 0.006 for 16x16 
FMIM and —0.502 ± 0.008 for 16 x 16 FFIM. To within error bars, the slopes seem to have an universal behavior 

max{Ai/fc} cx \og{fk)-^'^. (12) 

as predicted by Zhou and Bhatt . Our results suggest that the values of the slope is generic to the Wang-Landau 
algorithm and does not depend on system sizes and models. Certainly many more simulations on different models are 
needed to confirm the universality of the slope. 

IV. EFFECTS OF MODIFICATION FACTORS 

We also looked at how the Wang-Landau algorithm performs with different sequences of modification factors. In 
the extreme case, we studied the effects of taking the limit of / = 1 only after a few iterations. Assuming n iterations 
were performed with large modification factors, and on the {n + l)th iteration, the modification factor is set to 1. 
The background for implementation is as follows: Eq. Q uses the Boltzmann weights {B{E,T) = exp{—E/kBT)) 
and the resulting energy distribution is, 

PiE) = g*{E)B{E,T)/Z = g* {E)cM- E / k^T) / Z (13) 

where g*{E) is the true density of states. In general other weights can be used in summing the partition function. If 
one chooses B{E,T) = l/g„(£'), then the probability distribution of E for the {n + l)th iteration P„_|_i(i?) will be 
given by, 

P„+i{E)^g*{E)/g^{E)Z (14) 

where gn{E) is the density of states calculated by the nth iteration. Z is is an undetermined normalization constant. 
The true density of states can then be estimated by the accumulated histogram of the (n -\- l)th iteration. 

gn+i{E) = Hn+i{E)gn{E) X constant (15) 
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FIG. 3: Comparison of accuracy of the Wang-Landau method with two modification sequences, 
with the sequence fk+i = y/Jk (empty circles) and with the sequence where the limiting value of / : 
(filled squares). Lattice size is 32 x 32 with energy range E/Ne G [—1.55, —1.35]. 



Data points were generated 
= 1 is used after 14 iterations 



This is analogous to the iteration process employed in Lee's entropic sampling but we used the fast diffusion 
of the Wang-Landau algorithm in the early stage. Fig. |2| compares the accuracy of the Wang-Landau method with 
two different modification sequences for the 32 x 32 ferromagnetic Ising model. The energy range was E/Ne G 
[—1.55, —1.35] where Ne — 1024 is the total number of lattice sites. The vertical axis is the error of density of states 
defined by, 

771 ^ 

E 

where g* (E) is the exact density of states calculated from a MATHEMATICA program provided by Beale . gn+i{E) 
is the calculated density of states and m is the total number of energy levels in the summation over this energy range. 
We plot A for different sequences of modification factors. Empty circles were generated with modification factors 
/fc+i — \fjk (with /o = exp(l)) and stopping when the condition {H^^^~ H^^^) / {H^^^+ H^^^) < 0.1 is satisfied. Where 
i?,-„ax and i?,„i„ are the maximum and minimum histogram counts respectively. Filled squares were obtained with a 
sequence of modification factors where the limiting value of / = 1 was used after 14 iterations. The arrow indicates 
the location which the modification factor was set to 1. We measure the errors (filled squares) at several Monte 
Carlo steps per site after we set f — 1- Error bars were obtained by averaging over several independent simulations. 
Accuracy increases rapidly immediately after setting / = 1, but in the long run, the limiting case becomes only about 
twice as accurate as the other sequence. 



97i+i{E) 



9*{E) 



(16) 



V. CONCLUSION 



We derived an expression for the error term of the Wang-Landau algorithm. With this, we showed that the 
fluctuation of the accumulated histogram AHk plays a central role in the accuracy of the Wang-Landau method. We 
have also proposed that stopping each iteration as soon as AHk becomes saturated would be optimal. The dependence 
of the saturation values on the modification factor was also investigated and it was found that for the ferromagnetic 
Ising model and fully frustrated Ising model, maxjAiJ^} cx log(/fc)~^/^. With this equation, the saturation values of 
AHk for small modification factors can be predicted from values obtained with larger modification factors. Perhaps, 
a more efficient algorithm can be developed. We also studied the effects of using different sequences of modification 
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factors (refinement), in which we presented the Hmiting case where the modification factor is set to 1 after 14 iterations. 
There are significant improvements of efficiency for short simulations and improvements become less for longer runs. 
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